Emergence of spatiotemporal chaos driven by far-field breakup of spiral waves in the 

plankton ecological systems 
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Alexander B. Medvinsky et al [A. B. Medvinsky, I. A. Tikhonova, R. R. Aliev, B.-L. Li, Z.-S. Lin, 
and H. Malchow, Phys. Rev. E 64, 021915 (2001)] and Marcus R. Garvie et al [M. R. Garvie and G 
Trenchea, SIAM J. Gontrol. Optim. 46, 775-791 (2007)] shown that the minimal spatially extended 
reaction-diffusion model of phytoplankton-zooplankton can exhibit both regular, chaotic behavior, 
and spatiotemporal patterns in a patchy environment. Based on that, the spatial plankton model 
is furtherly investigated by means of computer simulations and theoretical analysis in the present 
paper when its parameters would be expected in the case of mixed Turing-Hopf bifurcation region. 
Our results show that the spiral waves exist in that region and the spatiotemporal chaos emerge, 
which arise from the far-field breakup of the spiral waves over large ranges of diffusion coefficients 
of phytoplankton and zooplankton. Moreover, the spatiotemporal chaos arising from the far- field 
breakup of spiral waves does not gradually involve the whole space within that region. Our results 
are confirmed by means of computation spectra and nonlinear bifurcation of wave trains. Finally, 
we give some explanations about the spatially structured patterns from the community level. 
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I. INTRODUCTION 



There is a growing interest in the spatial pattern dy- 
namics of ecological systems i, i, i, i, [1,0, i, B m 
[TTI , IT2I , However, many mechanisms of the spatio- 

temporal variability of natural plankton populations are 
not known yet. Pronounced physical patterns like ther- 
moclines, upwelling, fronts and eddies often set the frame 
for the biological process. Measurements of the underwa- 
ter light field are made with state-of-the-art instruments 
and used to calculate concentrations of phytoplankton 
biomass (as chlorophyll) as well as other forms of organic 
matter. Very high diffusion of the marine environment 
would prevent the formation of any stable patch spatial 
distribution with much longer life-time than the typical 
time of biodynamics. Meanwhile, in addition to very 
changeable transient spatial patterns, there also exist 
other spatial patterns in marine environment, much more 
stable spatial structure associated with ocean fronts, spa- 
tiotemporal chaos [10, [O, [ll]? cyclonic rings, and so 
called meddles [l5|. In fact, it is significant to create 
the biological basis for understanding spatial patterns 
of plankton [16]. For instance, the impact of space on 
the persistence of enriched ecological systems was proved 
in laboratory experiments [l^- Recently, it has been 
shown both in laboratory experiments [18] and theoreti- 
cally [13, E, H, mi that the existence of a spatial struc- 
ture makes a predator-prey system less prone to extinc- 



* Corresponding author; Electronic address: j inzhn@263.net| 



tion. This is due to the temporal variations of the density 
of different sub-populations can become asynchronous 
and the events of local extinction can be compensated 
due to re-colonization from other sites in the space |22|]. 
During a long period of time, all the spiral waves have 
been widely observed in diverse physical, chemical, and 
biological systems [23|, 0, [25I. [2^ . However, a quite lim- 
ited number of documents [ul, [l3, [13, [H, [1^ concern 
the spiral wave pattern and its breakup in the ecological 
systems. 

The investigation of transition from regular patterns 
to spatiotemporally chaotic dynamics in spatially ex- 
tended systems remains a challenge in nonlinear sci- 
ence [l3, 123I , [30! , [31] . In a nonlinear ecology system, the 
two most commonly seen patterns are spiral waves and 
turbulence (spatio-temporal chaos) for the level of the 
community [32] . It has been recently shown that sponta- 
neous spatiatemoporal pattern formation is an instrinsic 
property of a predator-prey system [llj, [13, [H, HI, |35|, [s^ 
and spatiotemporal structures play an important role in 
ecological systems. For example, spatially induced speci- 
ation prevents the extinction of the predator-prey mod- 
els [ll|, |l2|, [33 • So far, plankton patchiness has been ob- 
served on a wide range of spatial temporal scales [38, 39] . 
There exist various, often heuristic explanations of the 
spatial patterns phenomenon for these systems. It should 
be noted that, although conclusive evidence of ecological 
chaos is still to be found, there is a growing number of 
indications of chaos in real ecosystems [40|, [4l|, [42|, [43] . 

Recently developed models show that spatial self- 
structuring in multispecies systems can meet both cri- 
teria and provide a rich substrate for community-level 
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section and a major transition in evolution. In present 
paper, the scenario in the spatiahy extended plankton 
ecological system is observed by means of the numeri- 
cal simulation. The system has been demonstrated to 
exhibit regular or chaostic, depending on the initial con- 
ditions and the parameter values [lO, • We find that 
the far-field breakup of the spiral wave leads to complex 
spatiotemporal chaos (or a turbulentlike state) in the spa- 
tially extended plankton model ([1]). Our results show 
that regular spiral wave pattern shifts into spatiotempo- 
ral chaos pattern by modulating the diffusion coefficients 
of the species. 



II. MODEL 

In this paper we study the spatially extended nutrient- 
phytoplankton-zooplankton-fish reaction-diffusion sys- 
tem. Following Scheffer's minimal approach 0], which 
was originally formulated as a system of ordinary diff- 
ential equation (ODEs) and later developed models [TO, 
[Til l29|, l45|, |46|, as a further investigation, we study a 
two-variable phytoplankton and zooplankton model on 
the level of the community to describe pattern formation 
with the diffusion. The dimensionless model is written 
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where the parameters are r, a, 6, m, n, dp, (i^, and 
/ which refer to work in Refs. |Tl|. The explana- 
tion of model ([1]) relates to the nutrient-phytoplankton- 
zooplankton-fish ecological system [see Refs. |lo|, HH, 13 
for details]. The local dynamics are given by 
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From the earlier results [45i] about non-spatial system 
of model ([1]) by means of numerical bifurcation analysis 
show that the bifurcation and bistability can be found in 
the system ([1]) when the parameters are varied within a 
realistic range. For the fixed parameters (see the caption 
of Fig. [T]and[2|), we can see that the / controls the dis- 
tance from Hopf bifurcation. For larger /, there exists 
only one stable steady state. As / is decreased further, 
the homogeneous steady state undergoes a saddle node 
bifurcation (SN), that is fsN = 0.658. In this case, a 
stable and an unstable steady state become existence. 
Moreover, the bistability will emerge when the parame- 
ter / lies the interval fsN>f>fc = 0.445 (this value is 



more than the Hopf onset, fn = 0.3397). There are three 
steady states: with these kinetics A and C are linearly 
stable while B is unstable. Outside this interval, the sys- 
tem has unique nontrivial equilibrium. Recent stud- 
ies [11LI29I shown that the systems ^ can well-develop 
the spiral waves in the oscillation regime, but where the 
authors only consider the special case, i.e., dp = dh. A 
few important issue have not yet been properly addressed 
such as the spatial pattern if dp dh. 

Here we report the result that emergence of spatiotem- 
poral chaos due to breakup in the system under the 
dh 7^ dp case. We may now use the / and diffusion 
ratio, v = dh/dp^ as control parameters to evaluate 
the region for the spiral wave. Turing instability in 
reaction-diffusion can be recast in terms of matrix sta- 
bility [13, E^. Such with the help of Maple software 
assistance algebra computing, we obtain the parameters 
space (/, u) bifurcation diagrams of the spiral waves as 
showing Fig. O in which two lines are plotted, Hopf line 
(solid) and Turing lines (dotted) respectively. In domain 
I, located above all three bifurcation lines, the homo- 
geneous steady states is the only stable solution of the 
system. Domain II are regions of homogeneous oscilla- 
tion in two dimensional spaces [10]. In domain HI, both 
Hopf and Turing instabilities occur, (i.e., mixed Turing- 
Hopf modes arise), in which the system generally pro- 
duces the phase waves. Our results show that the system 
has spiral wave in this regions. One can see that a Hopf 
bifurcation can occur at the steady when the parameter 
/ passes through a critical values fn while the diffusion 
coefficients dp = dh = and the bifurcation periodic so- 
lutions are stable. From our analysis (see Fig. [2]), one 
could also see that the diffusion can induce Turing type 
instability for the spatial homogeneous stable periodic 
solutions and the spatially extended model ^ exhibit 
spatio-temporal chaos patterns. These spatial pattern 
formation arise from interaction between Hopf and Tur- 
ing modes, and their subharmonics near hte codimension- 
two Hopf- Turing bifucation point. Special, it is interest- 
ing that spiral wave and travelling wave will appear when 
the parameters correspond to the Turing-Hopf bifurca- 
tion region HI in the spatially extended model ([1]), i.e., 
the Turing instability and Hopf bifurcation occur simul- 
taneously. 



III. NUMERICAL RESULTS 

The simulation is done in a two-dimensional (2D) 
Cartesian coordinate system with a grid size of 600 x 600. 
The fourth order Runger-Kutta integrating method is 
applied with a time step At = 0.005 time unit and a 
space step Ax = Ay = 0.20 length unit. The results 
remain the same when the reaction-diffusion equations 
were solved numerically in one and two spatial dimen- 
sions using a finite-difference approximation for the spa- 
tial derivatives and an explicit Euler method for the time 
integration. Neumann (zero-flux) boundary conditions 
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FIG. 1: The sketch map for the bistabihty and the Hopf bi- 
furcation in the system (|2]) with r = 5.0, a = 5.0, b = 5.0, 
m = 0.6, and n = 0.4. The black curve is the gi{p,h). The 
colored curves are g2{p, h) with different values of /. The red 
curve: / = 0.3; the blue: / = 0.445; the green: / = 0.5; and 
the cyan: / = 0.658. 
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FIG. 2: The sketch map of parameter space (/, z/) bifurcation 
diagrams for the spatially extended system ([1]) with r = 5.0, 
a 5.0, h 5.0, m 0.6, 0.05, and n 0.4. 



were emmployed in our simulation. The diffusion terms 
in Eqs. (p!a|) and (pJp]) often describe the spatial mixing 
of species due to self-motion of the organism. The typi- 
cal diffusion coefficient of plankton patterns is about 
0.05, based on the parameters estimatie of Refs [50, 51] 
using the relationship between turbulent diffusion and 
the scale of the space in the sea. In the previous stud- 
ies [10, 11, 29, 45, 46], the authors provided a valueable 
insight into the role of spatial pattern for the system ([T]) 
if dp = dh' From the biological meaning, the diffusion 
coefficients should satisfy dh > dp. However, in nature 
waters it is turbulent diffusion that is supposed to domi- 
nate plankton mixing [52], when dh < dp is allowed. The 
other reason for choosing such parameter is that it is well- 



known new patterns, such as Turing patterns, can emerge 
in reaction-diffusion systems in which there is an imbal- 
ance between the diffusion coefficients dp and dh 0, [H^ . 
Therefore, we set u = dh/dp^ and investigated whether a 
spiral wave would break up into complex spatiotemporal 
chaos when the diffusion ratio was varied. Throughout 
this paper, we fix (ip = 0.05 and dh is a control parameter. 

In the following, we will show that the dynamic behav- 
ior of the spiral wave qualitatively change as the control 
parameter dh increases from zero, i.e., the diffusion ra- 
tio V increases from zero, to more than one. For large 
V {u > 1), the outwardly rotating spiral wave is com- 
pletely stable everywhere, and fills in the space when the 
proper parameters are chosen, as shown in Fig. [31(A). Fig- 
ure m A) shows a series of snapshots of a well-developed 
single spiral wave formed spontaneously for the variable 
p in system ([T]). The spiral is initiated on a 600 x 600 grid 
by the cross-field protocol (the initial distribution chosen 
in the form of allocated "constant-gradient" perturbation 
of the co-existence steady state) and zero boundary con- 
ditions are employed for simulations in the two dimen- 
sions. From Fig. [31(A) we can see that the well-developed 
spiral waves are formed firstly by the evolution. Inside 
the domain, new waves emerge, but are evolved by the 
spiral wave growing from the center. The spiral wave 
can steadily grow and finally prevail over the whole do- 
main (a movie illustrating the dynamical evolution for 
this case [54] [partly movie-1^ movie-2^ and movieS 
for dh = 0.2]). Fig. [S^B) shows that the spiral wave 
first break up far away from the core center and even- 
tually relatively large spiral fragments are surrounded 
by a 'turbulent' bath remain. The size of the surviv- 
ing part of the spiral does not shrink when dh is further 
decreasing until finally dh equals to 0, which is different 
from phenomenon that is observed previous in the two- 
dimensional space Belousov-Zhabotinsky and FitzHugn- 
Nagumo oscillatory system [30|, [Sll, [55|, |5y, 113 , in which 
the breakup gradually invaded the stable region near the 
core center, and finally the spiral wave broke up in the 
whole medium. Figure [31(C) is the time sequences (ar- 
bitrary units) of the variables p and h at an arbitrary 
spatial point within the spiral wave region, from which 
we can see that the spiral waves are caused by the ac- 
cepted as "phase waves" with substantially group veloc- 
ity, phase velocity and sinusoidal oscillation rather than 
the relaxational oscillation with large amplitude. This 
breakup scenario is similar to the breakup of rotating 
spiral waves observed in numerical simulation in chemi- 
cal systems 30^ |3l|, HH, HH, 113 , and experiments in BZ 
systems [58|, ^|, which shows that spiral wave breakup 
in these systems was related to the Eckhaus instability 
and more important, the absolute instability. 

The corresponding trajectories of the spiral core and 
the spiral arm (far away from the core center) at y = 300 
are shown in Fig. [H respectively. From Fig. [H we can 
see that the spiral core is not completely fixed, but oscil- 
lates with a large amplitude. However, as dh decreases 
to a critical value, an unstable modulation develops in 
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FIG. 3: Well developed spiral waves and some properties of 
them. The figures show simulations of the system ([1]) with 
r = 5, a = 5, 6 = 5, m = 0.6, n = 0.4, dp = 0.05, and 
/ = 0.3. (A) Well developed spiral waves shown at subsequent 
snapshot in time, = 0.2. (B) Far- field breakup of the spiral 
waves shown at subsequent snapshot in time, dh = 0.002. 
The white (black) areas correspond to maximum (minimum) 
values of p [Additional movie format available from Ref. [54]]. 
(C) Oscillations of the variable p and h at an arbitrary spatial 
point within the regular spiral wave region for both scenarios. 
Each figure is ran the long time until it spatial patterns are 
unchange. 



regions which is far away from the spiral core (cf. the 
middle column of the Fig. [4j). These oscillations eventu- 
ally grow large enough to cause the spiral arm far away 
from the core to breakup into complex multiple spiral 
waves, while the core region remains stable (the corre- 
sponding movie can be viewed in the online supplemen- 
tal in Ref. [sS] [partly movie -1 and movie -2^ and for 
dh = 0.02]). Figures [3l[B) and [D^B) show the dynamic 



behavior for dh = 0.02, i.e., u = 0.4. The regular tra- 
jectories far away from the core are now the same as the 
region of the spatial chaos (cf. the middle column of the 
Fig. [4|). It is shown that an decrease in the diffusion ra- 
tio u which leads to population oscillations of increasing 
amplitude (cf. the left column of the Fig. [4]). In the 
tradition explain that the minimum value of the popula- 
tion density decreases and population extinction becomes 
more probable due to stochastic environmental perturba- 
tions. However, from the spatial evolution of system ^ 
(see Fig. [3]), the temporal variations of the density of 
different sub-population can become asynchronous and 
the events of local extinction can be compensated due to 
re-colonization (or diffusion) from other sites. 




FIG. 4: The corresponding trajectories (from left to right) 
for locations (300,300), (250,300), and (50,300) respectively. 
The parameters in (A), and (B) were the same as these in 
Fig. [3j A) and (B), respectively. 

Furthermore, it is well known that the basic arguments 
in spiral stability analysis can be carried out by reducing 
the system to one dimensional space [sol, [111, [HH, HH, [131 • 
Here we show some essential properties of the spiral 
breakup resulting from the numerical simulation. In the 
next section we will give the theoretical computation by 
using the eigenvalue spectra. In this model, it is worth 
noting that we do not neglect the oscillation of the dy- 
namics in the core as shown in Fig. 2] due to the system 
exhibiting spatial periodic wave trains when the model 
is simulated in one-dimensional space. Breakup occurs 
first far away from the core (the source of waves). The 
spiral wave breaks towards the core until it gets to some 
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constant distance and then the surviving part of the spi- 
ral wave stays stable. These minimal stable wavelengths 
are called Xmin- So the one-parameter family may be 
described by a dispersion curve X{dh) (see Fig. [5]). The 
minimal stable wavelength Xmin of the spiral wave are 
shown in Fig. [5] coming from the simulation in two di- 
mensional space. The results of Fig. [5] can be interpreted 
as follows: the minimal stable wavelengths decrease with 
respect to the decrease of dh but eventually stay at a 
relative constant value, which is that the stable spiral 
waves are always existing for a larger region values of d^. 
Space-time plots at different times are shown in Fig. [6] 
for two different dh^ i.e., different which display the 
time evolution of the spiral wave along the cross section 
in the two-dimensional images of Fig. [SjA) and (B). As 
shown in Fig. E^A) and (B) for dh = 0.2 and dh = 0.02 
respectively, the waves far away from the core display 
unstable modulated perturbation due to convective in- 
stability [30|, I31I5 l55|, |56|, |57|, but this perturbation is 
gradually advected to the left and right sides, and finally 
disappears. The instability manifests itself to produce 
the wave train breakup several waves from the far-field, 
as shown in Figs.[6](B). 
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FIG. 5: Dependence of the wavelength Xmin on the parameter 
dh for the system ([!]) with r — 5.0, a — 5.0, b — 5.0, m — 0.6, 
dp — 0.05, and n — 0.4. Note the log scale for du- 




FIG. 6: Space-time plots of variable p for different time and 
dh- The parameters in (A), and (B) are the same as those in 
Fig. [3jA) and (B), respectively. 

spectrum is also required to compute, since it determined 
only by the far-field wave trains of the spiral. The lin- 
ear stability spectrum consists of point eigenvalues and 
the essential spectrum that is a continuous spectrum for 
spiral waves. 

For sake of simplicity, the Eqs. ([Taj) and (pJp]) can been 
written as following 

dp 

— - = dpV'^p + gi{p, h), (3a) 



IV. SPECTRA AND NONLINEAR 
BIFURCATION OF THE SPIRAL WAVE 

In this section, we concentrate on the linear stabil- 
ity analysis of spiral wave by using the spectrum the- 
ory [Hi, 0, EH, [62|, m . From the results in Kefs. 0, El 
we know that the absolute spectrum must be computed 
numerically for any given reaction-diffusion systems. In 
practice, such computations only require discretization 
in one-dimensional space and compare with computing 
eigenvalues of the full stability problem on a large do- 
main due to the spiral wave exhibitting traveling waves 
in the plane (see Fig. [6] about the space-time graphes). 
For spiral waves on the unbounded plane, the essential 



— =dhV^h^g2{v,h). (3b) 

Suppose that (p*, /i*) are a solutions and refer to them 
as steady spirals of Eq. ([3]) that rotate rigidly with a 
constant angular velocity cj, and that are asymptotically 
periodic along rays in the plane. In a coratating coordi- 
nate frame, usingthe standardized analysis method for 
the spiral waves [62|, [H] , the Eq. (|3]) is given by 

% = dp^leP + + ^ib*, ^*), (4a) 



— = dh^loh + + ^2(P*, /^*), (4b) 
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where (p, 0) denote polar coordinates, spirals waves are 
relative equilibria, then the statianry solutions p*{p^O) 
and 9) both are 27r-periodic functions with = cp — 
cot. In Eqs. (|4a|) and (|4bp the operator q denotes dpp-\- 



A. Computation of spiral spectra 

Next, we commpute the leading part of its linear stabil- 
ity spectrum for the system (j4]). Consider the linearized 
evolution equation in the rotating frame, the eigenvalue 
problem of Eqs. (|4a|) and (|4b|) associated with the planar 
spiral solutions 0) and 0) are given by 



dp 



dpVloP + + ^1 (P*, + 9i{p". h'')h = Xp, (5a) 



dh 

dhVleh + ^^+ gl{p\h*)p + g^(p*,h*)h = \h, (5b) 

where • • • , ^2 denote the derivatives of the nonlin- 
ear functions and gi{p^h) = r{l — p) — rp 



l+6p 



j^, gliP, h) = -T^,. 9l{p. h) = ^^- and 
52^P, h) = ^^-m-^ + We shall ignore 

isolated eigenvalues that belong to the point spectrum, 
instabilities caused by point eigenvalues lead to mean- 
der ingor drifting waves, or to an unstable tip motionin 
in excitable media and oscillation media [H, IH, IH, . 
This phenomenon is not shown in the present paper. In- 
stead, we focus on the continuous spectrum that is re- 
sponsible for the spiral wave breakup in the far field (see 
Fig. [31(b)). By the results in Ref. [62], it turns out that 
the boundary of the continuous spectrum depends only 
on the limiting equation for p ^ 00. Thus, we have that 
A is the boundary of the continuous spectrum if, and only 
if the limiting equation 

dpK,pP + 9l(jP*,h*)p + 5? (p*, h*)h = Xp, (6a) 



have solutions p(p, S) and /i(p, S) for (p, S) G x [0, 27r], 
which are bounded but does not decay as p ^ 00. Since 
spiral waves are rotating waves in the plane, the wave 
train solutions have the form as u(t^ x, y) = u{p^ cp — cot) 
for an appropriate wave numbers k and temporal fre- 
quency cj, where we assume that u is 27r-periodic in 
its argument so that u{(_) = u{(_ -\- 27r) for all ^ and 
u = (p^h)^. Spiral waves converge to wave trains 
u{p^(p — cot) Uyjtikp + (/? — cot) for p ^ 00, which 
are corresponding to asymptotically Archimedean in the 
two-dimensional space. Assume that A: 7^ and 7^ 0, 
and in this case, we can pass from the theoretical frame p 



to the comoving frame ^ = kp-\-(p — cot (^ G 7^) in which 
the eigenvalue equation (|6]) becomes 

dpk'^Wl^p^cop^^g{{u^t{0)P^9i{uwt{0)h = Ap, (7a) 



dhk^Vl^h + coh^ + g^{u^t{0)P + 92{uwt{0)h = Xh{7h) 

Indeed, any nontrivial solution u{^) = {p{£,) ^ h{^))^ cor- 
responding to the linearization eigenvalue problem ([7j) 
give a solution U{p^-) of the eigenvalue problem for the 
temporal period map of (|3|) in the corotating frame via 

U{p, •) = e^'u{kp - cot), U{p, T) = e^^u{kp - 2it). 

(8) 

We write the equations ([7j) as the first-order systems 

dp 
dp 
dh 
dp 
dp_ 
dp 
dhi 
dp 



Pi, 



hi, 



- = k '^dp^[pp-u;pi- gl{u^t{0)p- 9i{uwt{0)h]^ 



k '^d^ ^ [ph - cjhi - g2{uwt{0)P - 92{uwt{0)h] , 

(9) 



in the radial variable p. Then the spatial eigenvalues or 
spatial Floquet exponents are deternined as the roots of 
the Wronskian 



/ 



A{\,k) :-- 



1 



V 



1 



where k e U. The function U{p, •) = e^^e'^Puo{kp 






(10) 
cot) 

satisfies the equation (|3]) when the spatial and temporal 
exponents \k and A satisfy the complex dispersion rela- 
tion det(^(A,A:) - iA:) = for A G C. We cah the \k 
in spectrum of A(A, k) as spatial eigenvalues or spatial 
Floquet exponents. 

The stability of the spiral waves state (p*, /i*) on the 
plane is determined by the essential spectrum given by 

^ess = {A G C; det(^(A, k) -\k) ={) for some k G U}. 

(11) 

Now, we compute the continuous spectrum with the 
equation ([9]) that are parameterized by the wave num- 
ber k. For each A, there are infinitely many stable and 
unstable spatial eigenvalues. We plot A in the complex 
plane associated spatial spectrum, see Fig. [71 By the ex- 
plaination of Sandstede et al [HO], one would know that 
if the real part of the essentail spectra is positive, then 
the associated eigenmodes grow exponentially toward the 
boundary, i.e., they correspond to a far-field instability. 
Note that we find the essentail spectra are not sensitive 
to temporal frequency, co. 
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velocity is given by 



FIG. 7: The essentail spectra of wave trains are obtained by 
using the algorithms outlined in Refs. [goI . IgJ]. The param- 
eters of (A) and (B) are corresponding to the values used in 
the simulations of Fig. [3jA) and (B). 



B. Existence and properties of wave trains 

Suppose that a reaction-diffusion system on the one- 
dimensional space such that the variables equal to a 
homogeneous stationary solution. If the homogeneous 
steady-state destabilizes, then its linearization accommo- 
dates waves of the form e^^^^"'^^) for certain values k and 
CO. Typically, near the transition to instability, small spa- 
tially periodic travelling waves arise for any wave number 
close to /cc, which is the critical wavenumber. Their wave 
speed is approximately equal to where cOc is corre- 
sponding to kc. In present paper, we focus exclusively on 
the situation where ujc = and kc 0. The bifurcation 
with uJc = and kc is known as the Turing bifur- 
cation, and the bifurcating spatially periodic steady pat- 
terns are often referred to as Turing patterns. Another 
class of moved patterns will appear when the instabilities 
modulated by Hopf- Turing bifurcation, which is resem- 
ble a travelling waves. Moreover, the common feature 
of the spiral waves in one-dimensional space mentioned 
above is the presence of wave trains which are spatially 
periodic travelling waves of the form p^t{kx — cot; k) and 
hwt{kx — cot;k)^ where Pwt{(t^]k) and hyjt{4>]k) are 27r- 
periodic about (j). Typically, the spatial wavenumber k 
and the temporal frequency co are related via the non- 
linear dispersion relation co = co{k) so that the phase 



uo 



(12) 



A second quantity related to the nonlinear dispersion 
relation is the group velocity, Cg = of the wave 
train which also play a central role in the spiral waves. 
The group velocity Cg gives the speed of propagation of 
small localized wave-package perturbations of the wave 
train [g^]- Here, we are only concerned the existence of 
travelling wave solution. In fact, the spiral waves move 
at a constant speed outward from the core (see Fig. [6]), 
so that they have the mathematical form p(x, t) = P{z), 
and h{x,t) = H{z) where z = x — Cpt. Substituting these 
solution forms into Eq. (|3|) gives the ODEs 



(PP dP 



(13a) 



(13b) 



Here, we investigate numerically the existence, speed 
and wavelength of travelling wave patterns. Our ap- 
proach is to use the bifurcation package Matcont 2.4 [6S,] 
to study the pattern ODEs (p!3|) . To do this, the most 
natural bifurcation parameters are the wave speed Cp and 
/, but they give no information about the stability of 
travelling wave as solutions of the model PDEs (|3]) . 

Our starting point is the homogeneous steady state of 
Eq. ([13]) with in the domain HI of Fig. [2 The typical bi- 
furcation diagrams are illustrated in Fig. [HI which shows 
that steady spatially peroidic travelling waves exist for 
the larger values of the speed Cp, but it is unstable for 
small values of Cp. The changes in stability occur via 
Hopf bifurcation, from which a branch of periodic orbits 
emanate. Note that here we use the terms "stable" and 
"unstable" as referring to the ODEs system ([T3|) rather 
than the model PDEs. Fig. [8]^B) illustrates the max- 
imun stable wavelength against the bifucation parame- 
ter, speed Cp, and the small amplitudes have very long 
wavelength. It is known that Cp = ^, hence the tavelling 
wave solution exist when the Cp 7^ 0, i.e., /c 7^ 0,cc; 7^ 0. 
Using Matcont 2.4 package, it is possible to track the lo- 
cus of the Hopf bifurcation points and the Limit point 
(fold) bifurcation in a parameter plane, and a typical ex- 
ample of this for the Cp-f and Cp-dh plane are illustrated 
in Fig. [9l The travelling wave solutions exist for values 
of Cp and / lying in left of Hopf bifurcation locus (see 
Fig. [njA)). The same structure about the Cp-dh plane 
is shown in Fig. [H^B). These reuslts confirm our previ- 
ous analysis coming from the algebra computation (see 
Fig. [2]) and the numerical results (see Fig. [6]). 



V. CONCLUSIONS AND DISCUSSION 

We have investigated a spatially extended plank- 
ton ecological system within two-dimensional space and 
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FIG. 8: Typical bifurcation diagrams for the pattern 
ODEs ([13)). (A) The spatiahy periodic travehing waves of 
system JS]) is existence. The changes in stabihty occur via 
Hopf bifurcation, from which a branch of periodic orbits em- 
anate. Thus unstable travelling waves appear. (B) Maxi- 
mum stable wavelength along the bifurcation parameter Cp, 
i.e., /c / 0,0; / 0. The parameter values in (A) and (B) are 
the same as Fig. EJA). 
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found that its spatial patterns exhibit spiral waves dy- 
namics and spatial chaos patterns. Specially, the sce- 
nario of the spatiotemporal chaos patterns arising from 
the far-field breakup is observed. Our research is based 
on numerical analysis of a kinematic mimicking the dif- 
fusion in the dynamics of marine organisms, coupled to a 
two component plankton model on the level of the com- 
munity. By increasing (decreasing) the diffusion ratio of 
the two variables, the spiral arm first broke up into a 
turbulence-like state far away from the core center, but 
which do not invade the whole space. From the previous 
studies in the Belousov-Zhabotinsky reaction, we know 
the reason causing this phenomenon can be illuminated 
theoretically by the M. Bar and L. Brusch [sj], as 
well as by using the spectrum theory that poses by B. 
Sandstede, A. Scheel et al [Hi, [13, El, [H . The far-field 
breakup can be verified in field observation and is useful 
to understand the population dynamics of oceanic ecolog- 
ical systems. Such as that under certain conditions the 
interplay between wake (or ocean) structures and bio- 
logical growth leads to plankton blooms inside mesoscale 
hydrodynamic vortices that act as incubators of primary 
production. From Fig. [3] and corresponding the movies, 
we see that spatial peridic bloom appear in the phyto- 
plankton populations, and the details of spatial evolution 
of the distribution of the phytoplankton population dur- 
ing one bloom cycle, respectively. 

In Ref. [70], the authors study the optimal control of 
the model ([T]) from the spatiotemporal chaos to spiral 
waves by the parameters for fish predation treated as a 
multiplicative control variable. Spatial order emerges in a 
range of spatial models of multispecies interactions. Un- 
surprisingly, spatial models of multispecies systems often 



8.5 1 1.5 2 

(B) ^P^^^'S 

FIG. 9: An illustration of the variations in parameter space of 
the pattern ODEs ()13p . We plot the loci of Hopf bifurcation 
points. (A) f — Cp planes; (B) dh — Cp planes. The parameter 
values in (A) and (B) are the same as Fig. [3{A). 

manifests very different behaviors from their mean-field 
counterparts. Two important general features of spatial 
models of multispecies systems are that they allow the 
possibility of global persistence in spite of local extinc- 
tions and so are usually more stable than their mean-field 
equivalents, and have a tendency to self-organzie spa- 
tially or regular spatiotemporal patterns [tqI iTlj . The 
spatial structures produces nonrandom spatial patterns 
such as spiral waves and spatiotemporal chaos at scales 
much larger than the scale of interaction among individ- 
uals level. These structures are not explicitly coded but 
emerge from local interaction among individuals and lo- 
cal diffusion. 

As we know that plankton plays an important role in 
the marine ecosystem and the climate, because of their 
participation in the global carbon and nitrogen cycle at 
the base of the food chain [72[. From the review a 
recently developed ecosystem model incorporates differ- 
ent phytoplankton functional groups and their competi- 
tion for light and multiple nutrients. Simulations of these 
models at specific sites to explore future scenarios sug- 
gest that global environmental change, including global- 
warming-induced changes, will alter phytoplankton com- 
munity structure and hence alter global biogeochemical 
cycles [74]. The coupling of spatial ecosystem model to 
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global climate raises again a series of open questions on 
the complexity of model and relevant spatial scales. So 
the study of spatial model with large-scale is more impor- 
tant in the ecological system. Basing on numerical simu- 
lation on the spatial model, we can draft that the oceanic 
ecological systems show permanent spiral waves and spa- 
tiotemporal chaos in large-scale over a range of parame- 
ter values dh^ which indicates that periodically sustained 
plankton blooms in the local area. As with all areas of 
evolutionary biology, theoretical development advances 
more quickly than does empiraical evidence. The most 
powerful empirical approach is to conduct experiments 
in which the spatial pattern can be measured directly, 
but this is difficulties in the design. However, we can in- 
directly measured these phenomenona by the simulation 
and compared with the satellite pictures. For example, 
the spatiotemporal chaos patterns agree with the per- 



spective observatio n of the Fig. 3 in Ref. [73i]. Also, some 
satellite imageries [http: / / oceancol or.gsfc.nasa.gov| have 
displayed spiral patterns that represent the phytoplank- 
ton [the chlorophyll] biomass and thus demonstrated that 
plankton patterns in the ocean occur on much broader 
scales and therefore mechanisms thought diffusion should 
be considered. 
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